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PREFACE 

Work  reported  herein  was  performed  by  California  Research 
and  Technology,  Inc.,  (CRT)  during  the  period  June  1987  through 
October  1987  for  the  US  Army  Engineer  Waterways  Experiment 
Station  (WES)  under  Contract  No.  DACA39-84-D-0003  (Task  A-6,  Item 
0034).  It  was  sponsored  by  the  Office,  Chief  of  Engineers,  US 
Army,  as  a  part  of  Project  4A162784AT40 ,  Task  AO,  Work  Unit  026, 
"Blast-Induced  Soil  Liquefaction  Effects  on  Hardened  Structures 
and  Obstacle  Craters." 


Mr.  Paul  J.  Hassig  was  the  CRT  Principal  Investigator  for 
this  study  which  involved  development  of  a  one-dimensional 
implicit  numerical  technique  for  solving  the  dynamic  conservation 
equations  for  a  multiphase  (solid/fluid)  system,  including  the 
relative  flow  of  fluid  and  soil.  The  resulting  computer  code  is 
called  CRIME.  Mr.  Martin  Rosenblatt  and  Dr.  Douglas  Hatfield 
provided  invaluable  support  and  assistance  in  the  preliminary 
calculations  which  led  to  the  development  of  CRIME.  Mr.  Hassig 
prepared  this  report  with  the  help  of  Ms.  JoAnne  Clark,  who 
prepared  the  figures.  The  work  was  technically  supported  and 
monitored  by  Dr.  G.  Y.  Baladi  of  the  Geomechanics  Division, 
Structures  Laboratory  (SL),  WES.  Dr.  Y.  Marvin  Ito  was  the  CRT 
Project  Manager  for  Contract  No.  DACA39-84-D-0003 ;  Dr.  J.  G. 
Jackson,  Jr.,  was  the  WES  Contracting  Officer's  Representative. 
Mr.  Bryant  Mather  is  Chief,  SL. 

COL  Dwayne  G.  Lee,  CE,  is  Commander  and  Director  of  WES. 

Dr.  Robert  W.  Whalin  is  WES  Technical  Director. 
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CONVERSION  FACTORS,  NON-SI  TO 
SI  (METRIC)  UNITS  OF  MEASUREMENT 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted 
to  SI  (metric)  units  as  follows: 


Multiol 


To  Obtain 


feet 

0.3048 

metres 

feet  per  second  (ft/sec) 

0.3048 

metres  per 

second 

g's  (standard  free  fall) 

9.806650 

metres  per 

second  squared 

inches 

25.4 

millimetres 

kips  (force) 

4.448222 

kilonewtons 

kips  (force)  per  square  inch 

6.894757 

megapascals 

pounds  (mass) 

0.4535924 

kilograms 

AN  IMPLICIT  FINITE  DIFFERENCE  FORMULATION  FOR 
TREATING  MULTIPHASE  FLOW  IN  WET  POROUS  SOILS 


CHAPTER  1 

INTRODUCTION 

The  effects  of  multiphase  flow  on  the  cratering  and  ground 
motions  in  response  to  surface  bursts  (conventional  and  nuclear) 
in  layered  saturated  soils  is  not  yet  fully  understood.  The 
typical  state-of-the-art  cratering  calculations  are  just 
beginning  to  address  the  multiphase  aspects  with  regards  to  the 
constitutive  relationships  (material  models)  which  describe  the 
wet  soil  behavior.  Initial  attempts  describe  the  total  stress 
behavior  of  a  wet  soil  layer  and  the  effective  stress  behavior  of 
the  soil  lattice,  enabling  the  separate  determination  of  the 
pore-fluid  pressures  (Reference  1). 

The  M-DICE  code  was  developed  to  further  examine  the 
multiphase  issues  by  including  separate  models  of  the  pore-fluid 
pressures  and  soil  lattice  effective  stress  behavior.  More 
importantly,  the  relative  flow  of  the  pore-fluid  and  soil  lattice 
is  also  treated  by  calculating  the  separate  flow  of  each, 
including  their  mutual  interactions  due  to  drag.  With  this  fully 
multiphase  modeling,  the  M-DICE  code  was  used  to  calculate  the 
ground  motions  and  crater  formation  from  the  MISERS  BLUFF  II-l 
high-explosive  test  event  and  from  a  postulated  1  Mt  nuclear 
surface  burst  (References  2-5). 

Both  M-DICE  calculations  went  beyond  the  time  most  cratering 
codes  stop.  In  so  doing,  they  showed  the  existence  of 
significant  residual  pore  pressures  and  associated  relative 
velocities,  particularly  in  the  1  Mt  nuclear  calculation.  The 
subsequent  dissipation  and  redistribution  of  these  excess  pore 


$ 
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pressures  provide  a  potential  mechanism  for  significant  late  time 
ground  motions,  including  further  slumping  of  the  crater  slopes 
and  surface  settlements . 

Current  explicit  numerical  methodologies  are  not  practical 
for  calculating  to  very  late  times,  as  the  time  step  is  limited 
by  stress  wave  speeds .  Implicit  numerical  techniques  should 
remove  the  explicit  time  step  limitation,  allowing  for  much 
larger  time  steps  (factors  of  1000  or  more),  limited  perhaps  only 
by  the  particle  velocities.  Our  ultimate  objective  is  to  develop 
such  a  multiphase  implicit  formulation  in  two  dimensions  (2-D)  in 
order  to  efficiently  calculate  relative  flow  and  stress 
redistribution  in  wet  porous  soils.  The  focus  of  the  current 
study  is  the  development  of  a  suitable  technique  in  one  dimension 
( 1-D) .  The  resulting  finite-difference  code  is  called  CRIME 
(California  Research  Implicit  Multiphase  Eulerian) . 

The  implicit  multiphase  formulation  is  discussed  in  Chapter 
2  of  this  report.  Various  1-D  numerical  test  results  are 
presented  in  Chapter  3 .  A  summary  is  provided  in  Chapter  4  along 
with  conclusions  and  recommendations. 


CHAPTER  2 


METHODOLOGY 


2.1  GOVERNING  EQUATIONS 

The  following  equations  summarize  the  1-D  conservation 
equations  for  mass  and  momentum  for  a  solid-fluid  multiphase 
system  (energy  and  gravitational  considerations  are  ignored  for 
simplicity) .  The  equations  are  presented  in  differential  form. 
The  subscripts  f  and  s  refer  to  the  fluid  phase  and  solid  phase, 
respectively. 
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where  p  is  the  average  (or  bulk)  density, 

u  is  the  velocity, 

D  is  the  local  drag  force  coefficient, 

Pf  is  the  pore-fluid  pressure, 

PB  is  the  effective  mean  normal  stress  (pressure)  of  the 
solid  phase, 

Sx  is  the  effective  stress  deviator  of  the  solid  phase 
(negative  for  compression),  and 

u  is  the  local  area/volume  fraction  of  each  phase. 


The  total  stress  of  the  multiphase  material  is  modeled  in 
terms  of  the  effective  stress  of  the  solid  (soil  lattice)  and 
pore-fluid  pressure.  A  constitutive  model  is  used  to  describe 
the  effective  stress  behavior  of  the  solid  soil  lattice. 
Separate  hydrodynamic  equations  of  state  (EOS)  are  used  for  the 
water  and/or  air  to  describe  the  pore-fluid  pressures.  In  the 
case  of  partially  saturated  soils,  a  pressure  equilibrium  model 
is  used  for  the  water/air  mixture  EOS  in  the  pores.  Thus,  a 
fully  3-phase  model  involving  water-air-solid  is  used. 

The  following  equations  summarize  the  multiphase  equations 
of  state  for  the  fluid  phase  (pore  pressure)  and  solid  phase 
(effective  stress): 


Kf  dP f  Kf  (  dP f  +  1  P f  d pB 

P £  dt  P f  l  dt  pag  uf  dt  . 
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3  9xk  .  pB  dt 

?r  9u3  1  f  1 

al  3x  '  3  l  9xk  J  . 


=  2G 


=  2G< 


l_3  +  2.  dP& 


3  pg  dt 


where 


P  s  dt 


is  the  velocity  divergence. 
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Pf  =  Pf/^f  is  the  local  pore-fluid  density, 

i/f  =  1  -  vB  =  1  -  ps/pB g  is  the  porosity,  and 

Pag  is  the  average  density  of  the  individual  soil  grains 
(assumed  constant). 


The  solid  phase  EOS  (or  material  model)  describes  the 
effective  stress  behavior  of  the  solid  soil  lattice.  In  general 
the  bulk  modulus,  Ks,  and  shear  modulus,  GB,  are  local  functions 
of  position  and  time,  and  depend  on  the  nature  of  the  transient 
loading/unloading  behavior  of  the  soil  material.  Some  of  the 
numerical  examples  presented  in  this  study  assume  constant  values 
for  these  moduli.  Also  used  in  this  study  is  the  cap-type 
constitutive  model  of  Baladi  and  Barnes  (Reference  1).  This 
model  relates  incremental  changes  in  strain  (both  plastic  and 
elastic)  to  incremental  changes  in  stress. 

The  fluid  bulk  modulus,  K* ,  is  also  considered  to  vary  with 
position  and  time.  For  a  dry  soil,  an  adiabatic  ideal  gas 
equation  of  state  is  used  to  describe  the  pore  air,  such  that 

Pao 

Kf  =  Ka  =  7  Pa  -  (8a) 

Pa 

The  subscript  o  denotes  initial  values.  For  a  saturated  soil,  a 
hydro-elastic  equation  of  state  is  used  to  describe  the  water, 
such  that 


Kw  Pw 


Pw-Pwo 


For  a  partially  saturated  soil,  a  mixture  equation  of  state 
assuming  pressure  equilibrium  (Pa=pw)  is  used,  such  that 

Hi  -  Hi  +  Hi 

Kf  Ka  Kw 

The  drag  force  is  primarily  a  function  of  the  relative 
velocity,  ur  =  uf  -  us,  which  results  from  the  differential 
gradients  in  pore-fluid  pressure  and  effective  stress.  The 
relative  flow  of  the  fluid  with  respect  to  the  solid  is  also 
influenced  by  the  physical  properties  of  each  material,  i.e., 
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solid  phase  permeability,  porosity  and  particle  size,  as  well  as 
the  viscosity  of  the  pore-fluid.  In  this  study  a  modified 
Darcy's  law  is  used,  relating  the  apparent  volumetric  fluid  flow 
rate,  i/fur,  to  the  fluid  pressure  gradient,  i.e., 


f*  [  } 

=  -  ~  ^fur 
k  v  4 


where  is  the  fluid  viscosity  [g/cm/s]  and  k  is  the  local  soil 
permeability  [cm2].*  If  we  assume  a  locally  steady  state  and 
spatially  uniform  fluid  momentum  flux  in  Equation  3,  substitution 
of  Equation  9  into  Equation  3  yields  the  following  local  drag 


coefficient : 


=  a 

k 


The  permeability  is  assumed  to  vary  as  a  function  of  the  local 
porosity  (Reference  6),  i.e., 


°  [ufn )  J  k°  (1+e  )  lej 


where  e  =  vt/vB  is  the  void  ratio.  Figure  2.1  shows  the  behavior 
of  the  permeability  as  a  function  of  porosity  for  various  initial 
porosities . 

2 . 2  APPROACH 

The  finite-difference  analogs  of  the  time-dependent 
governing  equations  (1)  to  (7)  describe  the  complete  multiphase 
behavior  at  each  grid  location.  In  an  explicit  formulation  the 
terms  of  the  right  hand  side  (r.h.s.)  are  known  "time  n" 
quantities  and  are  used  to  obtain  the  left  hand  side  (l.h.s.) 
"time  n+1"  quantities  for  a  given  time  step.  The  resulting  time 
n+1  densities  and  velocities  are  then  used  to  get  the  final  time 


*  Permeability  units:  1  Darcy  =  9.87  x  10-9  cm2 

=  9.66  x  10~4  cm/s  (for  water) 


n+1  stresses  (pressures  and  stress  deviators),  thus  completing 
the  integration  over  one  time  cycle. 

In  general,  a  numerical  finite  difference  formulation  of  a 
time-dependent  set  of  equations  is  considered  "stable"  if  errors 
(truncation/round-off  errors,  etc.)  do  not  amplify  during  the 
time  integration  calculation.  Stability  is  usually  assured  by 
imposing  restrictions  on  the  size  of  the  time  step,  i.e.,  they 
are  "conditionally  stable."  For  typical  explicit  formulations, 
the  Courant-Friedrichs-Lewy  (CFL)  stability  criterion  requires 
the  time  step  (At)  to  be  less  than  the  ratio  of  the  grid  spacing 
(Ax)  to  the  stress  wave  speed  (c),  i.e., 

AtCFL  *  Ax  /  c  (12) 

Note  that  for  typical  cratering  and  ground  motion  calculations, 
the  stress  wave  speeds  are  determined  by  the  local  stress  moduli. 
These  wave  speeds  are  usually  much  larger  during  the  early  times 
(characterized  by  megabar  pressures)  compared  to  the  later  times, 
thus  allowing  for  increasing  time  steps  during  the  course  of  the 
calculation . 

Implicit  finite  difference  formulations  have  unknown  time 
n+1  values  on  the  r.h.s.  of  the  time-dependent  equations  in 
addition  to  the  l.h.s.  Usually  one  must  solve  a  set  of 
simultaneous  equations  to  obtain  the  updated  time  n+1  quantities 
(Reference  7).  The  advantage  of  such  formulations  is  to  improve 
the  stability  of  the  calculation  in  a  numerical  sense,  thus 
allowing  for  time  steps  much  larger  than  those  derived  from  the 
CFL  criterion.  The  disadvantage  of  such  schemes  is  the 
additional  amount  of  computations  often  needed  each  time  cycle  to 
solve  the  set  of  simultaneous  equations.  A  practical  implicit 
formulation  yields  a  significant  net  savings  in  computational 
time  compared  to  an  explicit  formulation,  i.e.,  the  decrease  in 
the  number  of  integration  cycles  realized  by  increasing  the  time 
step  outweighs  the  increase  in  computations  per  time  step. 
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Our  minimum  program  objective  is  the  development  of  a 
"practical"  implicit  formulation,  i.e.,  allowing  time  steps  at 
least  1000  times  greater  than  allowed  by  an  explicit  formulation. 
The  scheme  should  be  "conditionally  stable, "  with  time  steps 
restricted  only  by  the  material  particle  velocities,  u,  rather 
than  the  stress  wave  speeds,  i.e., 

Atu  <  Ax  /  u  (13) 

Thus,  in  a  manner  analogous  to  the  CFL  wave  speed  restriction, 
the  allowable  time  step  increases  as  the  particle  velocities 
decrease.  For  late  time  ground  motions  in  saturated  soils, 
typical  stress  wave  speeds  are  at  least  1500  m/s,  while  particle 
velocities  might  be  on  the  order  of  1-2  cm/s.  A  comparison  of 
the  particle  velocity  stability  requirement  (Equation  13)  and  the 
CFL  condition  (Equation  12)  shows  an  increase  in  allowable  time 
step  by  a  factor  of  100,0001  This  is  well  above  our  stated 
objective  of  at  least  1000. 

In  general,  the  choice  of  time  step  and  grid  size  for  a 
given  calculation  is  governed  by  the  desired  resolution,  or 
accuracy  of  the  solution.  For  example,  in  order  to  resolve  the 
propagation  of  stress  waves  for  a  given  grid  size,  the  time  step 
should  be  chosen  to  at  least  satisfy  the  CFL  criterion.  Larger 
time  steps  will  tend  to  diffuse  the  wave  fronts  over  time, 
similar  to  the  diffusion  which  occurs  for  coarser  grids. 
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2.3  FINITE  DIFFERENCING 


The  implicit  finite-differencing  scheme  of  the  time- 
dependent  governing  equations  (1)  to  (7)  is  presented  in  this 
section.  The  following  sketch  summarizes  the  primary  location  of 
the  physical  variables  in  the  1-D  x-direction: 


The  mass  (or  density)  cells  are  centered  at  i  (lower  case),  and 
the  momentum  cells  centered  at  I  (upper  case). 

The  fluid  and  solid  mass  conservation  Equations  (1)  and  (2) 
are  written: 

wr-  w:-M;-^[Nr-(p.«.);:;]  <« 
wr-  w:  =  w:-<[wr-  m:;i  <« 


where 


-  -  f1  -  ®)  sT"  [  (piur)  t  I1 

5  -  -  s»  si-  [M*)"  -  ('-•“•j °.j  <j 

The  time  centering  of  Equations  (14)  and  (15)  is  controlled  by 
the  variable  6.  A  value  of  0  =  0  implies  a  purely  explicit 


(14) 

(15) 

• 
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m 

finite  differencing,  while  a  value  of  0  =  1  implies  purely 
implicit  finite  differencing. 


The  fluid  and  solid  momentum  Equations  (3)  and  (4)  are 
written: 

wni 
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and  pj  = 


_  (bx-LPi  +  Axi  +  1pi  +  1 ) 


is  the  momentum  cell  density. 


(Axi  +  Axi+1 ) 

The  time  centering  of  Equations  (18)  and  (19)  is  controlled  by 
the  variables  X  and  <p  in  a  manner  analogous  to  6 . 

The  pore-fluid  pressure  Equation  (5)  and  effective  stress 
Equations  (6)  and  (7)  are  written: 
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In  order  to  more  readily  solve  the  finite  difference 
equations  for  the  fluid  and  solid  phase  momentums  at  each  grid 
location,  the  time  n+l  drag  and  solid  phase  stress  deviator  terms 
are  approximated  as  follows: 


n+l 

DX 
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(28) 
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(29) 


Substituting  Equations  (14)  and  (15)  into  Equations  (22),  (23) 
and  (24)  eliminates  the  unknown  time  n+l  densities.  Substituting 
the  resulting  equations  along  with  Equations  (28)  and  (29)  into 
the  momentum  equations  (18)  and  (19)  eliminates  the  time  n+l 
pressures  and  stress  deviators,  yielding  the  following  equations: 


oint,  governed  b; 
.t.  The  two  unkm 
,o  the  unknown  va 
esulting  set  of  i 
rid  are  obtained 


Table  2.1.  Definition  of  coefficients  evaluated  at  time  n  in 
fluid  momentum  equation  (30) 
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Table  2.2.  Definition  of  coefficients  evaluated  at  time  n  in 
solid  momentum  equation  (31) 
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CHAPTER  3 


NUMERICAL  RESULTS 

Several  different  1-D  (planar  geometry)  calculations  were 
performed  to  demonstrate  the  ability  of  CRIME  to  implicitly  treat 
various  aspects  of  the  multiphase  physics.  The  results  of  these 
test  cases  are  presented  in  four  subsections  emphasizing  the 
following: 

a)  Effective  stress  deviators  for  solid  phase  constitutive 
models  (velocity  and/or  strain  dependent,  e.g.,  cap- 
type), 

b)  Relative  flow  and  late-time  stress  redistribution 
(consolidation)  in  a  fluid-filled  porous  soil, 

c)  Comparison  with  an  independent  calculation,  and 

d)  Layered  geology  including  dry,  partially  and  fully 
saturated  soils. 

3.1.  EFFECTIVE  STRESS  DEVIATORS 

Table  3.1  summarizes  the  1-D  calculations  which  demonstrate 
the  ability  of  CRIME  to  implicitly  treat  effective  stress 
deviators.  A  constant  loading  pressure  of  50  bars  is  applied  to 
one  end  of  a  porous  soil  lattice  50  m  long,  which  is  divided  into 
50  uniform  cells  (Ax  =  1  m) .  A  rigid  (no-flow)  boundary 
condition  is  imposed  at  the  other  end. 

The  cap-type  constitutive  model  given  in  Reference  1  is  used 
to  describe  the  effective  stress  behavior  of  the  solid  soil 
lattice.  Figure  3.1  shows  the  uniaxial  strain  (UX)  load-unload 
stress-strain  behavior  for  vertical  strains*  up  to  12%  and  Figure 
3.2  shows  the  corresponding  effective  UX  stress  path  and  the 
effective  triaxial  shear  (TX)  failure  surface.  The  initial 


*  Strains  are  engineering,  i.e.,  evol  =  (vD  -  v)/vQ  =  (p  -  pQ)/p 


loading  bulk  modulus,  Ks  =  1.0  kbar,  and  shear  modulus, 

Gs  =  0.75  kbar,  yield  a  seismic  velocity  of  340  m/s  for  an 
initial  soil  density  pBO  =  1.75  gm/cc .  The  (explicit)  CFL 
condition  would  require  that  the  time  step  be  less  than 
AtCFL  =  2.9  msec.  The  signal  from  the  initial  loading  will  reach 
the  rigid  boundary  by  0.15  sec. 

Four  calculations  were  performed,  with  time  steps  varying 
from  0.5  msec  (0.17  times  A tCFL)  to  0.5  sec  (170  times  AtCFL) . 

The  implicit  capability  of  CRIME  is  exercised  by  setting  the 
numerical  parameters  controlling  the  time  centering  of  the  finite 
difference  equations  0  =  <f>  =  X  =  1  for  all  four  cases.  The 
results  are  presented  in  the  form  of  stress  and  velocity  profiles 
at  various  times. 

Figure  3.3  shows  the  total  effective  axial  stress*  and 
velocity  profiles  for  Case  E220,  which  used  a  time  step 
At  =  0.5  msec.  A  constant  50  bar  stress  wave  is  seen  to 
propagate  at  ~340  m/s,  with  an  associated  particle  velocity  of 
8.4  m/s.  The  shock  front  is  spread  over  ~5  zones,  the  midpoint 
of  which  is  at  33  m  by  50  msec.  Figure  3.4  shows  tha ^  the 
effective  pressure  and  axial  stress  deviator  components  are 
31  bar  and  19  bar,  respectively. 

Figures  3.5  and  3.6  show  the  stress  and  velocity  profiles 
for  Case  E223,  which  used  a  time  step  of  At  =  5  msec.  This  is  a 
ten-fold  increase  over  the  previous  case,  and  is  1.7  times  AtCFL. 
Although  the  calculation  needs  only  ten  time  steps  to  reach 
50  msec,  the  time  step  is  too  small  to  resolve  a  sharp  shock 
front;  the  larger  time  step  results  in  an  increased  smearing  of 
the  shock  front.  Note,  however,  that  the  general  characteristics 
of  the  solution  are  maintained,  as  the  midpoint  of  the  shock 
front  has  reached  ~33  m  at  50  msec. 


*  Stresses/pressures  are  plotted  positive  for  compression 
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Figures  3.7  and  3.8  show  the  stress  and  velocity  profiles 
for  Case  E224,  which  used  a  time  step  of  At  =  50  msec.  This  is 
an  additional  ten-fold  increase  over  the  previous  case,  and  is  17 
times  AtCFL.  In  this  case,  the  initial  stress  wave  will 
propagate  through  17  zones  in  one  time  step.  After  five  time 
steps  (0.250  sec),  the  wave  will  have  already  reflected  off  the 
left  boundary  and  propagated  2/3  of  the  way  back  to  the  source. 
Although  the  larger  time  step  has  resulted  in  a  further  smearing 
of  the  shock  front,  the  general  characteristics  of  the  solution 
are  still  maintained  for  each  of  the  five  integration  cycles,  and 
the  solution  has  remained  stable  in  a  numerical  sense. 

An  additional  calculation  was  performed  with  five  time  steps 
of  At  =  0.5  sec  (Case  E225,  Figures  3.9  and  3.10).  By 
t  =  1.5  sec,  the  calculation  is  essentially  complete,  as  the  the 
entire  soil  lattice  has  reached  zero  velocity  and  50  bars  of 
stress . 

3.2.  RELATIVE  FLOW  AND  LATE  TIME  STRESS  REDISTRIBUTION 

Table  3.2  summarizes  the  1-D  calculations  which  demonstrate 
the  ability  of  CRIME  to  implicitly  treat  relative  flow  and  late¬ 
time  stress  redistribution  (consolidation)  in  a  saturated  porous 
soil.  A  50-m-long  water-filled  soil  lattice  is  divided  into  50 
uniform  cells  (Ax  =  1  m) ,  with  rigid  boundaries  at  each  end. 

From  x  =  0  to  25  m  the  initial  water  pressure  is  1  bar  and  the 
effective  stress  is  2  bars;  from  25  m  to  50  m  the  water  pressure 
is  2  bars  and  the  effective  stress  is  1  bar.  Thus  the  initial 
total  stress  load  of  the  system  is  3  bars  everywhere. 

A  simple  linear-elastic  equation  of  state  is  used  to 
describe  the  effective  stress  behavior  of  the  soil  lattice,  with 
bulk  modulus  Kg  =  5  kbar  and  shear  modulus  Gs  =  3.75  kbar.  A 
simple  hydro-elastic  equation  of  state  is  used  to  describe  the 
water  pressure  behavior,  with  bulk  modulus  Kw  =  20  kbar. 
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The  rate  at  which  the  individual  stress  imbalances  reach 
equilibrium  (consolidation)  is  determined  by  the  soil 
permeability.  Two  calculations  were  performed  to  examine  this 
effect:  Case  E129  assumed  an  initial  permeability  of  100  cm/s , 

while  Case  E130  used  0.01  cm/s.  Each  calculation  started  with  a 
time  step  At  =  0.1  msec,  and  increased  by  a  factor  of  10  each 
subsequent  integration  cycle.  (Note  that  AtCFL  =  0.5  msec  for 
this  material.)  The  numerical  parameters  controlling  the  time 
centering  of  the  finite  difference  equations  were  set  to 
6  =  <p  =  X  =  1  for  each  case.  The  results  are  presented  in  the 
form  of  stress  and  velocity  profiles  at  various  times. 

Figure  3.11  shows  the  fluid  pressure  profiles  at  each  time 
(top  figure)  and  the  effective  vertical  stress  profiles  (lower 
figure)  for  Case  E129  (k  =  100  cm/s);  Figure  3.12  shows  the 
associated  fluid  and  solid  phase  velocity  profiles.  During  the 
consolidation  process,  peak  velocities  of  -7.4  cm/s  for  the  water 
and  +2.2  cm/s  for  the  solid  phase  occur  after  three  integration 
cycles  (t  =  11.1  msec).  These  peak  values  occur  at  the  point 
marking  the  position  of  the  original  stress  discontinuities 
(x  =  25  m),  and  thus,  the  largest  differential  stress  gradients. 
After  five  integration  cycles  (t  =  1.11  sec),  the  calculation  is 
essentially  complete,  as  both  the  fluid  pressure  and  effective 
stress  have  been  uniformly  redistributed  to  a  value  of  1.5  bar, 
and  the  particle  velocities  are  zero. 

When  the  permeability  is  reduced  by  a  factor  of  1000,  the 
consolidation  process  takes  longer,  as  demonstrated  by  Case  E130 
(k  =  0.01  cm/s).  Figures  3.13  and  3.14  are  similar  to  Figures 
3.11  and  3.12  for  Case  E129,  showing  the  first  five  integration 
cycles  to  t  =  1.11  sec.  Peak  particle  velocities  are  nearly  30 
times  smaller  and  the  consolidation  process  has  only  propagated 
to  ~10  m  on  either  side  of  x  =  25  m.  Figures  3.15  and  3.16  show 
the  final  three  integration  cycles  to  t  =  1111  sec,  at  which  time 
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the  calculation  is  considered  complete.  Note  that  the  time  for 
completion  of  the  consolidation  process  has  increased  by  a  factor 
of  ~1000,  consistent  with  the  decrease  in  soil  permeability. 

3.3.  COMPARISON  CALCULATION 

As  a  partial  validation  of  the  CRIME  code  a  calculation  was 
performed  which  can  be  compared  directly  to  that  performed  by 
Prevost  (Reference  8).  Table  3.3  summarizes  the  initial  and 
boundary  conditions  of  the  problem,  which  is  essentially  one  of 
1-D  wave  propagation  in  a  fluid-saturated  porous  medium. 

The  calculation  is  initiated  with  a  constant  10  bar  loading 
pressure  applied  at  the  surface.  10  m  of  dry  soil  (the  pores  are 
void)  overly  a  fluid-saturated  soil  lattice.  The  top  70  m  are 
divided  into  70  uniform  cells  (Ay  =  1  m) ;  subsequent  cells 
increase  by  10%  downward  to  over  1100  m  depth.  The  rigid  bottom 
boundary  is  chosen  far  enough  away  as  to  not  influence  the 
calculation  over  the  integration  time  frame  (1.7  sec).  At  the 
10  m  water  table  depth  a  partially  reflected  stress  wave  will 
propagate  back  to  the  surface,  where  computational  "absorbing 
dampers"  will  allow  it  to  pass  without  further  reflections. 

The  permeability  of  the  the  soil  lattice  is  assumed  to  be 
k  =  50  cm/s,  which  corresponds  to  a  drag  coefficient 
D  =  5  g/cc/s.  The  effective  stress  behavior  is  described  by  a 
simple  linear-elastic  equation  of  state,  with  bulk  modulus 
Ka  =  33.3  bar  and  shear  modulus  Ga  =  50  bar.  The  pore-fluid 
pressure  is  described  by  a  simple  hydro-elastic  EOS,  with  bulk 
modulus  Kf  =  5  kbar,  and  is  much  stiffer  than  the  soil  lattice. 

Both  CRIME  Case  E444  and  the  calculation  of  Prevost  used  a 
constant  time  step  At  =  8.5  msec.  Pore-fluid  pressures  and 
effective  stresses  are  compared  in  Figures  3.17  and  3.18  which 
show  time  histories  at  various  depths  in  the  saturated  region. 


Figures  3.19  and  3.20  show  comparisons  of  the  fluid  and  solid 
phase  particle  velocities.  The  results  are  quite  similar.*  The 
more  diffusive  nature  of  Case  E444  is  due  to  implementing  the 
fully  implicit  finite  differencing  of  CRIME  by  setting  the 
numerical  parameters  6  =  <p  =  X  =  1.  Prevost  implemented  a  "split 
operator  method",  whereby  the  solid  phase  equations  were  treated 
explicitly,  while  the  fluid  was  treated  implicitly.  The  exact 
duplication  of  Prevost' s  results  was  not  the  focus  of  this  study, 
and  thus  this  comparison  is  felt  to  be  adequate. 

3.4.  LAYERED  GEOLOGY  WITH  VARYING  SATURATION 

A  primary  objective  in  the  development  of  the  CRIME  code  was 
the  capability  of  treating  realistic  geologies.  Table  3.4 
summarizes  a  1-D  calculation  of  the  loading  and  subsequent 
consolidation  of  a  layered  geology  with  varying  saturation. 

The  site  profile  is  an  idealized  representation  of  the  top 
three  layers  of  the  MISERS  BLUFF  II-l  site  (Reference  1).  The 
initial  porosity  for  each  layer  is  assumed  to  be  27.4%:  Material 
A  consists  of  dry  sand  with  air-filled  porosity  ( 0  to  5  m  depth), 
Material  3  is  partially  saturated  (S=90%  water-filled  porosity,  5 
to  13  m  depth),  and  Material  C  is  fully  saturated  (13  to  34  m 
depth).  The  permeability  is  assumed  to  be  k  =  0.01  cm/s. 

The  effective  stress  behavior  described  by  the  cap-type  model 
of  Baladi  and  Barnes  (Reference  1)  is  assumed  the  same  for  each 
material.  Their  total  stress  behavior  differs  according  to  the 
degree  of  saturation,  as  demonstrated  in  Figure  3.21  which  shows 
the  UX  stress-strain  response  of  Materials  A,  B  and  C  in  terms  of 
total  stress.  Note  that  the  difference  between  the 


*  Note  that  the  sign  convention  of  Prevost  is  followed,  such  that 
compressive  stresses  are  negative,  compressive  pressures  are 
positive,  and  positive  velocities  are  downward. 
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Material  C  and  Material  A  stresses  represents  the  pore-water 
pressure  contribution  in  Material  C. 

The  computational  domain  consisted  of  1-m  cells  in  each 
layer.  A  constant  50  bar  loading  pressure  is  applied  to  the 
surface,  and  the  bottom  boundary  is  assumed  rigid.  The  initial 
pore  pressures  in  each  material  is  one  atmosphere  (1.0132  bar). 
Material  C  has  the  largest  seismic  velocity  (1900  m/s),  and  as 
such,  the  (explicit)  CFL  condition  would  require  that  the  time 
step  be  less  than  AtCFL  =  0.5  msec. 

Two  calculations  were  performed.  Case  E605  used  a  constant 
time  step  of  At  =  AtCFL  =  0.5  msec  to  a  maximum  simulation  time 
of  0.5  sec.  Case  E607  initially  used  a  time  step  of  At  =  5.0 
msec  to  a  simulation  time  of  0.6  sec.  Subsequent  integration 
runs  used  larger  time  steps,  with  an  eventual  maximum  of 
At  =  100  sec  from  750  to  1250  sec.  The  results  from  Case  E605 
are  presented  below,  followed  by  Case  E607. 

Figure  3.22  shows  the  total  vertical  stress,  and  Figure  3.23 
shows  the  corresponding  pore-fluid  pressure  and  effective 
vertical  stress  at  times  of  10,  20,  30,  40  and  50  msec  for  Case 
E605  (At  =  0,5  msec).  The  stress  wave  generated  by  the  50  bar 
overpressure  loading  propagates  downward,  arriving  at  the 
Material  A/B  interface  (5-m  depth)  after  10  msec.  The  total 
stress  behavior  between  these  two  materials  is  very  similar  at 
these  relatively  low  stress  levels  (Figure  3.21)  such  that  there 
is  very  little  reflection  occurring  at  this  interface.  However, 
the  partitioning  of  the  total  stress  loading  is  different. 

Whereas  virtually  all  of  the  50-bar  load  is  taken  by  the  soil 
lattice  in  Material  A,  the  water-air  mixture  in  the  pores  of 
Material  B  accounts  for  ~4  bars  of  the  total  stress. 

The  Material  B/C  interface  is  reached  after  30  msec.  The 
differing  total  stress  behavior  (impedance)  of  these  materials 
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causes  a  partially  reflected  stress  of  ~125  bar,  which  propagates 
in  both  directions.  Most  of  the  total  stress  transmitted  down 
through  Material  C  is  now  taken  by  the  pore-water  in  the 
saturated  soil.  By  50  msec  the  bottom  boundary  has  been  reached, 
generating  even  larger  reflected  stresses. 

Figure  3.24  shows  the  corresponding  velocities  of  the  soil 
lattice  and  relative  velocity  of  the  pore-fluid  with  respect  to 
the  soil  lattice.  During  the  loading  phase,  downward  velocities 
reach  ~8  m/s  in  Materials  A  and  B,  and  ~3  m/s  in  Material  C. 

Peak  upward  relative  velocities  of  ~17  cm/s  occur  at  the  Material 
B/C  interface.  It  is  at  this  interface  that  the  largest 
differential  pore  pressure/effective  stress  gradient  occurs  after 
the  passage  of  the  loading  stress  wave. 

Figures  3.25  to  3.29  show  the  multiphase  stress  and  velocity 
time  histories  (to  t  =  0.5  sec)  at  various  depths  for  Case  E605. 
The  propagation  of  the  initial  stress  loading  signal  and 
subsequent  reflections  and  transmissions  at  the  interfaces  and 
boundaries  are  evident  in  the  total  vertical  stress  time 
histories  (Figure  3.25).  The  maximum  reflected  stress  at  the 
bottom  boundary  of  nearly  250  bars  occurs  at  ~60  msec.  The 
constant  overpressure  loading  of  50  bars  at  the  surface  acts  to 
damp  out  the  incoming  reflected  stress  waves,  and  the  entire 
column  oscillates  around  the  51-bar  stress  level,  which  will  be 
the  final  stress  state. 

The  effective  stress  and  pore  pressure  behavior  are  shown  in 
Figures  3.26  and  3.27.  As  expected,  these  time  histories 
oscillate  in  phase  with  the  total  stress.  The  underlying 
distribution  of  the  total  stress  load  to  the  effective  stress  and 
pore  pressure  in  each  material  is  essentially  unchanged  during 
the  initial  0.5  sec,  as  the  effective  stress  carries  most  of  the 
load  in  Material  A  and  the  pore-water  does  so  in  Material  C. 
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The  solid  particle  velocity  and  relative  fluid  velocity  are 
shown  in  Figures  3.28  and  3.29.  These  time  histories  also 
oscillate  in  phase  with  the  total  stress.  However,  note  that 
with  respect  to  the  (oscillating)  soil  lattice  there  are  upward 
pore-fluid  velocities,  with  maximum  values  of  ~17  cm/s  occurring 
at  the  original  Material  B/C  interface  at  the  13-m  depth.  These 
relative  velocities  are  indicative  of  the  long-term  consolidation 
process  which  will  occur  as  the  pore  pressures  in  Materials  B  and 
C  dissipate  towards  pressure  equilibrium. 

Case  E607  demonstrates  the  ability  of  the  CRIME  code  to 
calculate  the  late-time  consolidation  process.  A  constant  time 
step  of  At  =  5  msec  (10  times  AtCFL)  is  used  for  the  first 
0.6  sec  of  the  calculation,  the  results  of  which  can  be  compared 
directly  to  Case  E605.  Figures  3.30,  3.31,  and  3.32  show  the 
multiphase  stress  and  velocity  profiles  at  10,  20,  30,  40  and 
50  msec.  Although  there  is  a  loss  of  temporal  resolution,  the 
general  characteristics  of  the  solution  are  maintained.  In 
particular,  the  differential  stress  loading  in  each  material  and 
the  relative  fluid  velocity  at  the  Material  B/C  interface  are 
reproduced. 

Figures  3.33  to  3.37  show  the  multiphase  stress  and  velocity 
time  histories  to  t  =  0.5  sec.  When  comparing  the  total  vertical 
stress  to  Case  E605  it  is  evident  that  the  larger  time  step  of 
Case  E607  has  served  to  filter  of  damp  out  the  oscillations 
associated  with  the  propagation  back  and  forth  of  the  reflected 
stress  wave.  By  t  =  0.5  sec  the  solution  is  relatively 
quiescent,  achieving  stress  and  velocity  values  consistent  with 
the  mean  values  of  Case  E605. 

Subsequent  time  integrations  of  Case  E607  used  even  larger 
time  steps  (see  Table  3.4),  enabling  the  calculation  to 
efficiently  proceed  to  t  =  1250  sec.  The  total  CPU  time  used  is 


one  third  the  amount  used  for  Case  E605,  while  simulating  to  2500 
times  farther  in  time. 


Figures  3.38  to  3.41  show  the  multiphase  stress  and  velocity 
time  histories  to  t  =  1250  sec,  during  which  time  consolidation 
occurs.  As  the  large  pore  pressures  in  Materials  B  and  C 
dissipate  slowly  due  to  the  relative  upward  flow  of  water,  the 
soil  lattice  consolidates  and  the  effective  stresses  increase  to 
the  50-bar  loading  stress.  Figures  3.42  summarize  this  process 
in  the  form  of  pore  pressure  and  effective  stress  profiles  at 
late  times.  Figure  3.43  shows  the  displacements  of  selected 
fluid  and  solid  tracers  during  the  calculation.  There  is  a  net 
upward  displacement  of  the  pore- fluid  and  a  net  downward 
displacement  of  the  soil  lattice  as  it  compresses  to  carry  the 
50-bar  loading  stress. 
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Table  3.1.  Summary  of  1-D  calculations  emphasizing  solid  phase 
stress  deviators 


/I 

l< - 

/I 

Soil  lattice 

!  < - 

constant  loading 

/I 

|< - 

pressure,  50  bar 

/I 

w/  voids 

l< - 

/I 

1  < - 

0  50 

Position  (m) 


zoning 

Ax  =  1  m 

initial  soil  density 

p 8  =  1.75  g/cc 

soil  porosity 

l/f  =  27.4  % 

seismic  velocity 

c  =  340  m/s 

maximum  explicit 

time  step 

^CFL  K  2.9  msec 

Case 

Time  step 
( sec ) 

Simulation 
time  (sec) 

Integration 

cycles 

E220 

0.0005 

0.05 

100 

E223 

0.005 

0.05 

10 

E224 

0.05 

0.25 

5 
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Table  3.2. 


Summary  of  1-D  calculations  emphasizing  relative  flow 
and  stress  redistribution 


Water-filled  soil  lattice 


/I 

1 

1 

/I 

/I 

water 

pressure  =  1 

bar 

1 

1 

water 

pressure  =  2 

bar  1 

1 

/I 

/I 

solid 

stress  =  2 

bar 

1 

1 

solid 

stress  =  1 

bar  1 

1 

25 

Position  (m) 


zoning 

Ax  = 

1  m 

elastic  soil 

Ks  = 

5.0 

kbar 

Gs  = 

3.75 

i  kbar 

initial  soil  density 

PB  = 

1.75 

g/cc 

soil  porosity 

Uf  = 

27.4 

% 

multiphase  wave  speed 

c  = 

2027 

m/s 

maximum  explicit 

time  step  At 

CFL  = 

0 . 5  msec 

Case 

E129 

E130 

permeability,  k 

100 

cm/ s 

.01 

cm/ s 

Maximum  time  step  in  calculation* 

1 

sec 

1000 

sec 

Simulation  time 

1.11 

.  sec 

1111 

sec 

Number  of  integration  cycles 

5 

8 

The  initial  time  step  At0  =  0.1  msec. 

For  each  integration  cycle  At  increased  by  a  factor  of  10. 
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Table  3.3.  Summary  of  1-D  calculation  emphasizing  comparison 
with  Prevost  (1985)  calculation* 


CRIME  Case  E444 


soil  lattice 
fluid-filled 


soil  lattice  I  <- 
I  <- 

w/  voids  I  <- 


constant  loading 
pressure,  10  bar 
(w/  absorbing 
dampers ) 


1126 


10 

Depth  (m) 


zoning  Ay  =  1  m  to  70  m  depth,  increasing  by  10%  each 

cell  thereafter 


Initial  Material  Properties 
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Table  3.4.  Summary  of  1-D  calculations  emphasizing  layered 
geology  with  varying  saturation 


/I 

1 

1 

|  < - 

/I 

Material  C 

1 

Material  B 

1 

Material  A 

|  < - 

constant 

/I 

1 

1 

|  < - 

loading 

/I 

Saturated 

1 

Partially- 

1 

|  < - 

over- 

/I 

soil 

1 

saturated 

1 

Dry  soil 

|  < - 

pressure. 

/I 

S=100% 

1 

soil  S=90% 

1 

S  =  0% 

1  < - 

50  bar 

/I 

1 

1 

1  < - 

34 

13 

5 

0 

Depth  (m) 


Initial  Material  Properties 


Material 

wet 

density 

(g/cc) 

dry 

density 

(g/cc) 

water 

content 

S  (%) 

Air  void 
content 

va  (%) 

Seismic 
velocity 
c  (m/'s) 

permea¬ 
bility 
k  (cm/s) 

A 

1.75 

1.75 

0 

27.4 

340 

.01 

B 

2.00 

1.75 

90 

2.74 

340 

.01 

C 

2.02 

1.75 

100 

0 

1900 

.01 

o  Geologic  profile  represents  idealized  MISERS  BLUFF  II-l  site 

o  Initial  porosity  is  assumed  to  be  i/{  =  27.4%  for  each  layer 

o  Effective  stress  behavior  for  each  layer  is  described  by  the 
cap-type  model  of  Baladi  and  Barnes  (Reference  1) 
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Table  3.4  (continued) 

Calculational  Parameters 
zoning  (Ay)  1  m 

multiphase  wave  speed  (c)  1.9  k 


Maximum  explicit  time  step 
( CFL  condition,  AtCFL  <  Ay/c) 


1.9  km/s 

0.5  msec 


Time  step 

Integration 

Simulation 

Cumulative 

cycle 

time 

CPU  time* 

E605 

0.0005 

0  -  1000 

0.  - 

0.5 

25.0 

E607 

0.005 

0  -  120 

0.  - 

0.6 

3.5 

0.5 

120  -  220 

0.6  - 

50.6 

6.0 

5. 

220  -  260 

50.6  - 

250.6 

7.0 

50. 

260  -  270 

250.6  - 

750.6 

7.5 

100. 

270  -  275 

750.6  - 

1250.6 

7.8 

Computer  used  is  ELXSI  6400 


X*v" 
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Figure  3.3.  Effective  axial  stress  and  particle  velocity  versus 
position  at  various  times  for  Case  E220. 
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Figure  3.26.  Effective  vertical  stress  versus  time  (0  <  t  <  0.5  sec) 
at  various  depths  for  Case  E605. 
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Figure  3.27.  Pore-fluid  pressure  versus  time  (0  <  t.  <  0.5  sec) 
at  various  depths  for  Case  E605 . 
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Figure  3.43.  Displacements  of  selected  fluid  and  solid  tracers 
versus  time  for  Case  E607. 


CHAPTER  4 


SUMMARY ,  CONCLUSIONS,  AND  RECOMMENDATIONS 

As  part  of  a  WES  study  of  blast-induced  soil  liquefaction 
effects,  the  M-DICE  code  was  developed  to  provide  a  rational 
methodology  for  calculating  multiphase  flows  in  wet  porous  soils. 
The  code  includes  separate  treatments  of  the  pore-fluid  pressures 
and  soil  lattice  effective  stress  behavior,  which  thus  yields  the 
total  stress  behavior.  More  importantly,  the  relative  flow  of 
the  pore-fluid  and  soil  lattice  is  also  treated  by  calculating 
the  separate  flow  of  each,  including  their  mutual  interactions 
due  to  drag.  With  this  fully  multiphase  modeling,  the  M-DICE 
code  was  used  to  calculate  the  ground  motions  and  crater 
formation  from  the  MISERS  BLUFF  II-l  high-explosive  test  event 
and  from  a  postulated  1  Mt  nuclear  surface  burst  (Reference  2-5). 

Both  M-DICE  calculations  went  beyond  the  time  most  cratering 
codes  stop.  In  so  doing,  they  showed  the  existence  of 
significant  residual  pore  pressures  and  associated  relative 
velocities,  particularly  in  the  1  Mt  nuclear  calculation.  The 
subsequent  dissipation  and  redistribution  of  these  excess  pore 
pressures  provide  a  potential  mechanism  for  significant  late  time 
ground  motions,  including  further  slumping  of  the  crater  slopes 
and  surface  settlements.  Unfortunately,  because  M-DICE  (like  all 
current  numerical  cratering  and  ground  shock  codes)  is  formulated 
in  an  explicit  manner,  calculating  to  very  late  times  is  simply 
not  practical.  The  time  step  is  limited  by  a  stability  criterion 
which  is  a  function  of  stress  wave  speeds.  An  implicit  numerical 
technique,  however,  would  remove  the  explicit  time  step 
limitation,  allowing  for  much  larger  time  steps  (factors  of  1000 
or  more),  limited  by  the  desired  accuracy  and  perhaps  only  by  the 
particle  velocities. 


In  this  study  a  1-D  implicit  multiphase  finite  difference 
formulation  which  calculates  the  relative  flow  and  dynamic  stress 
behavior  in  wet  porous  soils  was  developed  and  incorporated  into 
a  computer  code  called  CRIME.  Various  1-D  test  cases  are 
presented  which  demonstrate  the  ability  of  CRIME  to  efficiently 
calculate  to  very  late  times  with  time  steps  much  larger  ( factors 
>  1000)  than  permitted  by  standard  explicit  techniques.  In 
particular,  the  loading  and  subsequent  consolidation  of  realistic 
layered  geology  and  varying  saturation  has  been  successfully 
simulated  as  follows:  A  constant  50-bar  overpressure  loading  is 
applied  to  the  surface  of  the  top  three  layers  (to  34-m  depth)  of 
the  MISERS  BLUFF  II-l  site.  The  early  time  (t  <  0.5  sec) 
distribution  of  the  load  is  determined  by  the  degree  of 
saturation.  Most  of  the  load  is  carried  by  the  soil  lattice  in 
the  dry  Material  A  (0  to  5m)  and  by  the  pore-water  in  the 
saturated  Material  C  (13  to  34  m) .  The  late  time  redistribution 
(consolidation)  of  the  differential  stress  loads  is  governed  by 
the  soil  permeability  (0.01  cm/s).  The  CRIME  calculation  of  this 
consolidation  process  to  over  20  minutes  is  summarized  in  Figure 
3.42;  implicit  time  steps  as  large  as  100  sec  (20,000  times 
explicit)  were  used. 

The  1-D  results  presented  in  this  report  show  that  the 
implicit  approach  we  have  taken  can  be  used  to  efficiently 
calculate  relative  flow  and  dynamic  stress  behavior  in  wet  porous 
soils.  Extending  the  finite  differencing  of  the  governing 
equations  to  2-D  would  be  relatively  straightforward,  i.e.,  there 
will  be  four  unknowns  governed  by  four  equations  at  each  grid 
point  in  the  2-D  formulation  instead  of  the  two  unknowns  and  two 
equations  in  the  1-D  formulation.  A  computationally  efficient 
algorithm  which  solves  the  resulting  set  of  coupled  algebraic 
equations  will  have  to  be  developed,  and  the  more  complicated 
configuration  of  various  air/water/rolid  interface  boundaries 
will  have  to  be  implicitly  resolved. 


We  recommend  that  the  implicit  CRIME  formulation  be  extended 
to  2-D,  and  that  its  efficiency  and  effectiveness  for  calculating 
blast-induced  cratering  and  ground  shock  phenomenology  to  very 
late  times  be  demonstrated. 
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